
#include <fstream>
using std::ifstream;
using std::ofstream;
#include <iostream>
#include <sstream>
using std::istringstream;
#include <string>
using std::string;
#include <vector>
using std::vector;

#include "TMV.h"
#include "TMV_Sym.h"
#include "dbg.h"
using std::endl;

#define DATECODE "0913"

#define NAME "lambda"
// Add NAME to the output names (lnlike.NAME.dat and predlike.NAME.debug)

#define WMAP3


std::ostream* dbgout=0;
bool XDEBUG=false;

const double PI = 3.141592653589793;
const double DEG = PI/180.;
const double ARCSEC = DEG/3600.;
const double ARCMIN = DEG/60.;

const double map_minr = 1.0*ARCMIN;
//const double map_minr = 5.0*ARCMIN;
const double map_maxr = 40.0*ARCMIN;
const double var_minr = 40.0*ARCMIN;
const double var_maxr = 100.0*ARCMIN;
const double singthresh = 0.01;
const int datafrac = 10; // take every nth line of the data

#ifdef DILLUT
const double DILLUTSQ = DILLUT*DILLUT;
#endif

void ReadData(tmv::Vector<double>*& r, tmv::Vector<double>*& mean,
    tmv::SymMatrix<double>*& cov, size_t& nmapout)
{
  dbg<<"Start ReadData\n";

  const size_t nfields = 12;
  const string datafile = "/home/mjarvis/plots/m2.out.";
  const char* fieldcodes[nfields] = { 
    "A", "G", "H", "J", "K", "L", "M", "N", "Q", "R", "T", "X" };

  vector<vector<double> > rmap(nfields);
  vector<vector<double> > map(nfields);
  //vector<vector<double> > sigmap(nfields);
  vector<vector<double> > rvar(nfields);
  vector<vector<double> > var(nfields);
  //vector<vector<double> > sigvar(nfields);

  for(size_t k=0;k<nfields;k++) {
    string name = datafile + DATECODE + '.' + fieldcodes[k];
    ifstream fin(name.c_str());
    if (!fin) myerror("opening data file:  ",name.c_str());

    dbg<<"Reading file: "<<name<<endl;
    string line;
    int count=0;
    while (getline(fin,line)) {
      if (line[0] == '#') continue;
      count++;
      if (count < datafrac) continue;
      count = 0;
      istringstream linein(line);
      double r,mapsq,mxsq,junk,gamsq;
      double sigmap,siggamsq;
      linein >>r>>mapsq>>sigmap>>mxsq>>junk>>junk>>junk>>junk>>junk>>gamsq>>siggamsq;
      r *= ARCMIN;

#ifdef BMODE
      mapsq += (BMODE)*mxsq;
#endif

#ifdef DILLUT
      mapsq *= DILLUTSQ;
      gamsq *= DILLUTSQ;
      //sigmap *= DILLUTSQ;
      //siggamsq *= DILLUTSQ;
#endif

      if (r > map_minr && r < map_maxr) {
	rmap[k].push_back(r);
	map[k].push_back(mapsq);
	//sigmap[k].push_back(sigmap);
      }
      if (r > var_minr && r < var_maxr) {
	rvar[k].push_back(r);
	var[k].push_back(gamsq);
	//sigvar[k].push_back(sigvar);
      }
    }
    dbg<<"sizes = "<<rmap[k].size()<<", "<<rvar[k].size()<<endl;
    if (k>0) {
      Assert(rmap[k].size() == rmap[0].size());
      Assert(rvar[k].size() == rvar[0].size());
    }
  }
  size_t nmap = rmap[0].size();
  size_t nvar = rvar[0].size();
  size_t ntot = nmap + nvar;
  dbg<<"nmap = "<<nmap<<", nvar = "<<nvar<<endl;

  r = new tmv::Vector<double>(ntot,0.);
  mean = new tmv::Vector<double>(ntot,0.);
  cov = new tmv::SymMatrix<double>(ntot,0.);

  for(size_t k=0;k<nfields;k++) {
    for(size_t i=0;i<nmap;i++) {
      (*r)[i] += rmap[k][i];
      (*mean)[i] += map[k][i];
      for(size_t j=0;j<=i;j++) {
	(*cov)(i,j) += map[k][i] * map[k][j];
      }
    }
    for(size_t i=0;i<nvar;i++) {
      (*r)[nmap+i] += rvar[k][i];
      (*mean)[nmap+i] += var[k][i];
      for(size_t j=0;j<=i;j++) {
	(*cov)(nmap+i,nmap+j) += var[k][i] * var[k][j];
      }
      for(size_t j=0;j<nmap;j++) {
	(*cov)(nmap+i,j) += var[k][i]*map[k][j];
      }
    }
  }

  dbg<<"mean = "<<*mean<<std::endl;
  dbg<<"cov(0..3,0..3) = "<<cov->subSymMatrix(0,3)<<std::endl;

  *r /= double(nfields);
  *mean /= double(nfields);
  *cov -= double(nfields)*(*mean ^ *mean);
  *cov /= double(nfields*(nfields-1));
  for(size_t i=0;i<ntot;i++) Assert((*cov)(i,i) >= 0.);

  dbg<<"cov(0..3,0..3) = "<<cov->subSymMatrix(0,3)<<std::endl;

#ifdef BMODE
  *cov /= sqrt(BMODE*BMODE+1);
#endif

  nmapout = nmap;

}

int main(int argc, char* argv[]) 
{
#ifdef NAME
  string debugname = string("predlike.") + NAME + ".debug";
#else
  string debugname = string("predlike.debug");
#endif
  dbgout = new ofstream(debugname.c_str());
  dbgout = &std::cout;

#ifdef DILLUT
  dbg<<"DILLUT = "<<DILLUT<<endl;
#endif
#ifdef BMODE
  dbg<<"BMODE = "<<BMODE<<endl;
#endif
  dbg<<"map range = "<<map_minr/ARCMIN<<" - "<<map_maxr/ARCMIN<<endl;
  dbg<<"var range = "<<var_minr/ARCMIN<<" - "<<var_maxr/ARCMIN<<endl;

  tmv::Vector<double>* r;
  tmv::Vector<double>* data;
  tmv::SymMatrix<double>* Sigma;
  size_t nmap;
  ReadData(r,data,Sigma,nmap);
  size_t ntot = r->size();

  Sigma->divideUsing(tmv::SV);
  Sigma->saveDiv();
  dbg<<"Before SetDiv cov(0..3,0..3) = "<<Sigma->subSymMatrix(0,3)<<std::endl;
  Sigma->setDiv();
  dbg<<"After SetDiv cov(0..3,0..3) = "<<Sigma->subSymMatrix(0,3)<<std::endl;
  Sigma->svd().thresh(singthresh,dbgout);

  ofstream fout("Jarvis_Map.dat");

  fout<<"# The data values are (in order): \n";
  for(size_t i=0;i<nmap;i++) {
    fout<<"# M_ap("<<(*r)[i]/ARCMIN<<")\n";
  }
  for(size_t i=nmap;i<data->size();i++) {
    fout<<"# Var_TH("<<(*r)[i]/ARCMIN<<")\n";
  }
  fout<<"#\n";
  fout<<"# Our definition of the aperture mass and shear variance are: \n#\n";
  fout<<"# M_ap(R) = (9 Om_m^2 H0^3 / 8 pi c^3) int( l^2 W_Map(lR) P_kappa(l) )\n";
  fout<<"# Var_TH(R) = (9 Om_m^2 H0^3 / 8 pi c^3) int( l^2 W_Var(lR) P_kappa(l) )\n";
  fout<<"# with W_Map(eta) = eta^4/4. * exp(-eta^2)\n";
  fout<<"# and  W_Var(eta) = 4 J_1(eta)^2 / eta^2\n#\n";
  fout<<"# The measurements of M_ap(R) and Var_TH are made by integrating over the \n";
  fout<<"# correlation function:\n#\n";
  fout<<"# M_ap(R) = int_x=0..2( x^2 dlogx [T+(x) xi+(xR) + T-(x) xi-(xR)]/2 )\n";
  fout<<"# Var_TH(R) = int_x=0..2( x^2 dlogx S+(x) xi+(xR) )\n";
  fout<<"# with T+(x) = (x^4 - 16 x^2 + 32)/128 exp(-x^2/4) \n";
  fout<<"#      T-(x) = x^4/128 exp(-x^2/4) \n";
  fout<<"#      S+(x) = 1/Pi (4 arccos(x/2) - x sqrt(4-x^2) )\n";
  fout<<"#      xi+(r) = <e1(r1) e1(r1+r) + e2(r1) e2(r1+r)>\n";
  fout<<"#      xi-(r) = <e1(r1) e1(r1+r) - e2(r1) e2(r1+r)>\n";
  fout<<"#      (e1,e2 for each galaxy is measured relative to line connecting them)\n";
  fout<<"#\n";
  fout<<"# The first line below is the number of M_ap values ("<<nmap<<"),\n";
  fout<<"# then the number of Var_TH values ("<<data->size()-nmap<<").\n";
  fout<<"# The next line is the R for each data point, given in arcminutes.\n";
  fout<<"# The next line is the "<<data->size()<<" data values.\n";
  fout<<"# Then there are "<<data->size()<<" lines for the (symmetric) covariance matrix.\n";
  fout<<"# Note: The covariance matrix is singular (ie. the data are oversampled)\n";
  fout<<"#       so you will want to use a singular value decomposition to properly\n";
  fout<<"#       use it for division, etc.\n";
  fout<<"#       FYI, I calculate the first few singular values as: \n";
  fout<<"#       (relative to the largest value)\n";
  fout<<"#       ";
  for(size_t i=0;i<8;i++) fout<<Sigma->svd().getS()(i)/Sigma->svd().getS()(0)<<" ";
  fout<<"\n";
  fout<<"#       I would recommend using at most 6 singular values for any\n";
  fout<<"#       chisq calculation using this matrix.\n";

  fout<<nmap<<"  "<<(data->size()-nmap)<<std::endl;
  for(size_t i=0;i<data->size();i++) fout<<(*r)[i]/ARCMIN<<" ";
  fout<<"\n";
  for(size_t i=0;i<data->size();i++) fout<<(*data)[i]<<" ";
  fout<<"\n";
  for(size_t i=0;i<data->size();i++) {
    for(size_t j=0;j<data->size();j++) 
      fout<<(*Sigma)(i,j)<<" ";
    fout<<"\n";
  }

  return 0;
}
